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Tentative evidence for excess GeV-scale gamma rays from the galactic center has been corrobo¬ 
rated by several groups, including the Fermi collaboration, on whose data the observation is based. 
Dark matter annihilation into standard model particles has been shown to give a good fit to the 
signal for a variety of final state particles, but generic models are inconsistent with constraints from 
direct detection. Models where the dark matter annihilates to mediators that subsequently decay 
are less constrained. We perform global fits of such models to recent data, allowing branching frac¬ 
tions to all possible fermionic final states to vary. The best fit models, including constraints from 
the AMS-02 experiment (and also antiproton ratio), require branching primarily to muons, with 
a ~ 10 — 20% admixture of b quarks, and no other species. This suggests models in which there 
are two scalar mediators that mix with the Higgs, and have masses consistent with such a decay 
pattern. The scalar that decays to A 1+ A t_ must therefore be lighter than 2 m T = 3.6 GeV. Such a 
small mass can cause Sommerfeld enhancement, which is useful to explain why the best-fit annihi¬ 
lation cross section is larger than the value needed for a thermal relic density. For light mediator 
masses (0.2 — 2) GeV, it can also naturally lead to elastic DM self-interactions at the right level for 
addressing discrepancies in small structure formation as predicted by collisionless cold dark matter. 


1. INTRODUCTION 

Indirect detection may offer the best hope of discov¬ 
ering the nature of dark matter beyond its generic grav¬ 
itational effects. Currently there is a strong hint from 
Fermi Large Area Telescope data of dark matter anni¬ 
hilation in the galactic center (GC), giving rise to an 
excess of gamma rays peaking at energies of several GeV 
m- The morphology is consistent with that expected 
from dark matter annihilations, and the required cross 
section is of the right order of magnitude for yielding 
the correct thermal relic density. Millisecond pulsars are 
the main conventional astrophysical explanation that has 
been proposed wm ( see however refs. [T21 ITT] for al¬ 
ternative explanations). Arguments against the pulsar 
explanation have been given in refs. [HUTTTI . While the 
debate continues, the dark matter explanation remains a 
possibility that is still being widely explored El-El- 

The simplest models have heavy s-channel media¬ 
tors leading to annihilations XX ff where / repre¬ 
sents standard model particles that lead to gamma rays 
through their decays (or possibly inverse Compton scat¬ 
tering). However even if the mediator does not couple 
directly to light quarks, its coupling to any charged par¬ 
ticle induces mixing with the photon at one loop, leading 
to dark matter scattering on nucleons. In many cases 
the induced coupling exceeds that allowed by direct de¬ 
tection searches, when the annihilation cross section for 
XX ff is large enough to explain the GC gamma-ray 
excess. Cosmic ray antiproton data also put pressure on 
direct annihilation to bb HUSH ( see however HQ)- 

To escape direct detection and antiproton constraints, 


one can consider the possibility of heavy mediators <fi that 
are produced on shell in the annihilations XX ~t </></>, and 
subsequently decay. Even relatively slow decays due to 
very small couplings of <j> to standard model particles 
can be consistent with the morphology of the GC signal, 
while making the models compatible with direct searches 

m-m- 

In the present study we reconsider light mediator mod¬ 
els, motivated by new data sets for the GC excess, one 
from ref. [B] (hereafter CCW) and the other presented 
by the Fermi collaboration itself [521 • These results al¬ 
low for heavier dark matter producing gamma rays up 
to higher energies ~ 100 GeV than the earlier determi¬ 
nations of ©111 which preferred lighter ~ 30 GeV dark 
matter, though the error bars are large enough for over¬ 
lap of the allowed parameters. In addition we take into 
account constraints on annihilations producing electrons, 
from the AMS-02 experiment [51] [55] . These turn out to 
strongly disfavor direct annihilations into electrons that 
would produce a bump-like feature in the electron spec¬ 
trum. 

We find that the CCW and Fermi data show a strong 
preference for annihilation into mediators that decay 
with branching fractions ~ 80 — 90% to muons and 
~ 10 — 20% to 6-quarks. From a theoretical model build¬ 
ing perspective this at first looks peculiar; for a single 
mediator to have such couplings is not suggested by any 
symmetry principle. On the other hand, singlet scalars 
that mix with the Higgs boson will decay preferentially 
into the heaviest possible standard model particle, since 
the couplings are just proportional to those of the Higgs. 
Therefore a more natural explanation is to have two me- 
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diators, one whose mass is between 2m^ and 2 m T , and 
the other with mass between 2 mb and 2 mt- (For a medi¬ 
ator with mass greater than 2 one must also consider 
decays into Higgs bosons.) By this combination of fit¬ 
ting to data and theoretical motivations, we are led to 
consider models with two scalar mediators with a hier¬ 
archy of masses. Remarkably, the best fit to the data 
puts the mediator masses in the ranges described above, 
which need not have been the case. This provides further 
experimental motivation for building models with several 
mediators, which at first might have seemed like a large 
theoretical leap. 

We begin our analysis with an agnostic view concern¬ 
ing theoretical models, assuming only a single mediator 
for simplicity, and allowing for arbitrary branching ra¬ 
tios into different lepton and quark flavors, whose spectra 
are provided by ref. [523 (hereafter PPPC). We describe 
the construction of the y 2 functions for fitting to the 
CCW and Fermi data, respectively, combined with those 
of AMS-02. The preference for final states consisting of 
an admixture of muons and b quarks is thereby estab¬ 
lished. We then focus on these final states, adding a sep¬ 
arate mediator coupled to each one, and refit to the data. 
In the following section, a simple model of two scalars 
mixing with the Higgs is presented. We subsequently 
check whether the preferred models are consistent with 


complementary constraints from direct detection, ther¬ 
mal relic density, cosmic microwave background, dwarf 
satellite observations, the antiproton flux ratio, and dark 
matter self-interactions. A simple two-mediator model is 
constructed that can fit the observations without exces¬ 
sive fine-tuning; Sommerfeld enhancement of the annihi¬ 
lation in the galaxy plays an important role. The model 
is shown to be consistent with LHC constraints on extra 
scalars mixing with the Higgs boson. 


2. PREDICTED SIGNAL 


The photon flux (GeV/cm 2 /s/sr) from annihilation of 
Majorana DM particles in a given region around the GC 
can be expressed as 
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where ( av ) is the averaged cross section for DM annihila¬ 
tion into two mediators, is the branching ratio 

for decays of the mediator (j> into // final states, and J 
is the averaged J factor for the region of interest (ROI), 
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For Dirac DM, one should replace the prefactor 1/2 —» 
1 /4 in eq. 0 , resulting in a cross section that is twice as 
large as that for Majorana DM, for a fixed observed flux. 
Unless otherwise stated we will assume Majorana DM in 
the following. The J factor depends upon the assumed 
shape of the DM density profile, which is commonly taken 
to be of the generalized NFW type, 


P(r) = p s 


jr/r s ) 7 
(1 + r/r s ) 3 * * * ~'i 


(3) 


Here we take r s = 20 kpc and p s such that p Q = 0.4 
GeV/crn 3 in the solar neighborhood, with Rq = 8.5 kpc. 
Ref. [i] finds that the morphology of the GC excess is 
best fit by taking 7 = 1.26, while ref. ’5] adopts 7 = 1 . 2 . 

dNl/dE 1 is the photon spectrum from a single DM 
annihilation into // final states, and comes from boosting 
the spectrum due to decay of the mediators <f> [52] : 

dN* _ 2 r E i x + dE\ dNl 

dE 7 (x+ - x-) J E ^ X _ E\ dE\ 

where dN/ /dE 1 is the photon spectrum from </> —» // in 
the rest frame of <f> and x± = m x /m ( f ) ± \J ( m x /m ^) * 2 — 1 . 


For most final states, a Monte Carlo generator is needed 
to predict the distribution of photons from hadronization 
and decays. We take the spectra dN$ /dE from PPPC 
(ref. [SB]), which is valid for mediator masses down to 10 
GeV. For lighter mediators, extrapolation of the PPPC 
results would be required, introducing inaccuracies into 
the predictions. However our subsequent fits will dictate 
the need only for lighter mediators decaying into muons, 
for which we use an analytic expression for dN['/ / dE 1 

[87] . given in appendix [B] The factor of 2 in Q accounts 
for the two mediators produced in the annihilation. 


3. DATA SETS 

We consider three data sets for the GC excess: those of 

CCW [S], Fermi [S3], and Daylan et al. [J. Detailed de¬ 

scriptions follow; the data are summarized in fig. [T] We 

further include the AMS-02 measurement of the spec¬ 
trum of cosmic ray electrons/positrons, and antiproton 
ratio data from the BESS, CAPRICE and PAMELA ex¬ 
periments. The antiproton data is not included in our 
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Figure 1: The three data sets for the GC gamma ray excess. 
The Fermi flux presented is the total flux from the 15° x 15° 
square around the GC; the other two fluxes are normalised 
accordingly with same DM profile. Error bars for CCW are 
taken from the diagonal components of their covariance ma¬ 
trix. Solid curves are the predictions of the best-fit models 
described in section [4] 


overall % 2 function; rather we will verify after fitting the 
data that our models are compatible with the antiproton 
constraints. 


3.1. CCW spectrum 


The CCW spectrum (ref. 0) is downloadable from 
[ 88 ) . where a covariance matrix S for computing the y 2 
for fits to the data is also provided. Then 
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where the sum is over 24 energy bins, and (E -1 )^- are the 
matrix elements of the inverse of the covariance matrix, 
obtained from columns 29-52 of the data file from 188]). 
This accounts for the correlations between the different 
energy bins. 

The ROI for CCW is a ±20° square around the GC; 
in galactic coordinates, where £ is the longitude and b is 
the latitude, the region is 


\t\ < 20 ° and 2 ° < |b| < 20 ° ( 6 ) 


Figure 2: Spectrum for GeV excess from Fermi data, ex¬ 
tracted from ref. [83] 

3.2. Fermi spectrum 

The Fermi collaboration has not officially released its 
data, but we have digitized it from the presentation in ref. 
[ 88 ] , and list the results in table |III| appendix [A] Fluxes 
are given in 20 energy bins, equally spaced in log 10 (i?) 
between 1 MeV and 89 GeV. We estimate the statistical 
errors from taking \Z~N for the number of total events in 
each bin and applying this to the part of the signal inter¬ 
preted as the excess. Ref. [S3] gives two characterizations 
of the spectrum of excess events, one which is presumed 
to be a power law in energy with an exponential cutoff, 
and the other being a separate fit to the excess in each 
bin. We adopt the latter for our analysis. 

In addition to the statistical errors, there is systematic 
uncertainty associated with assumptions about the tem¬ 
plates for background photons from pulsars and OB stars. 
We define the signal as the median between the upper 
and lower envelopes found from varying these templates, 
and the systematic error as the difference. This is added 
in quadrature with the statistical error to estimate the 
total uncertainty. The result is plotted in fig. [2] showing 
that errors are systematics-dominated at low energy, but 
mostly statistical at high energy. The \ 2 function is then 
defined in the usual way. 

The ROI for the Fermi data is a 15° x 15° square 
around the GC. Numerically integrating © we obtain 
J = 1.07 x 10 23 GeV 2 /cm 5 , again assuming generalized 
NFW parameter 7 = 1.2 and p Q = 0.4GeV/cm 3 . Here 
we do not average over the solid angle (giving J instead 
of J) because the Fermi data are reported as a total flux 
rather than an intensity flux. 


where the ± 2 ° region in latitude is masked out to remove 
the galactic disk. 

With these inputs and 7 = 1.2 for the generalized 
NFW distribution, the J factor in ([I]) is J = 2.062 x 
10 23 GeV 2 cm -5 for the CCW data. (This number agrees 
with [58].) 


3.3. Daylan et al. spectrum 

Although our primary purpose is to explore the impli¬ 
cations of the first two data sets which are more recent, 
for completeness we also apply our methodology to the 
GC excess spectrum determined in ref. [4] . It can be read 
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Figure 3: Electron spectra from dark matter annihilation into 
2e, 4e, 2/x and '2b, showing why annihilations directly into 
electrons are most strongly constrained by AMS-02 data. 

from fig. 5 of that paper, from which it is straightforward 
to define \ 1 2 f° r a given predicted flux. The spectrum 
shown there has been normalized to correspond to a J 
factor that is not averaged over any solid angle, but in¬ 
stead is defined along a line of sight 5° away from the 
GC. For 7 = 1.2, we thus find J = 9.09 x 10 23 GeV 2 /cm 5 
(again using p 0 = 0.4GeV/cm 3 * ). 

3.4. AMS positron/electron spectrum 

AMS-02 published their recent results on the measure¬ 
ment of the e + /e~ ratio and the separate electron and 
positron fluxes in [MJ , confirming a positron excess 
above the expectations from standard cosmic ray prop¬ 
agation scenarios that was previously seen by PAMELA 
and Fermi. The more precise measurements of AMS-02 
do not tell us whether the excess of the positron flux orig¬ 
inates from dark matter or pulsars, but the smooth spec¬ 
trum up to 0(100) GeV already puts strong constraints 
on dark matter models [TSUHSHUI- If light dark matter 
annihilates directly into e + e _ , even though the spectrum 
is initially a delta function of energy, after propagation it 
leads to bump-like feature in the observed spectra. The 
feature is more spread out in models with decays into 
mediators, \X ~► ( f > 4 > 4e, where the initial spectrum is 

box-like; nevertheless the final shape is still localized in 
energy and can be strongly constrained. The differences 
in shape of the electron spectra for relevant final states 
are illustrated in fig. [3] 

To simplify our analysis and to make the results more 
model-independent, we assume that only dark matter 
annihilating to e + e - or 2 e + 2 e^ can be constrained by 
current AMS-02 data, since these have harder and more 
localized spectra compared to other channels. The latter 
also generate positrons and electrons, but they come from 
three- or four-body final states, which make the spectrum 
much softer and more broad. It would be easy to absorb 


such nondistinctive contributions into the smooth back¬ 
ground . 1 Following m , we use polynomial functions to 
fit the logarithm of AMS spectra as the background. To 
obtain the electron/positron flux from light dark mat¬ 
ter, a cosmic ray propagation model is chosen with a 
relatively large magnetic field, in order to achieve con¬ 
servative bounds, since the energy loss is essentially due 
to the magnetic field and will soften the spectra. 

Before including the light dark matter contribution, 
the absolute % 2 of the background to fit against the AMS 
positron ratio data is y 2 = 35.06. This value does not de¬ 
pend on assumptions about the propagation model or so¬ 
lar modulations since we use fitting functions, rather than 
considering specific propagation models. However for the 
dark matter signal, the propagation and solar modulation 
uncertainties have to be considered. We marginalize over 
the range of the effective potential of solar modulation, 
[0,1 GV] and a relatively large magnetic field is chosen 
to propagate positrons and electrons from dark matter 
annihilation . 2 The 3cr exclusion curves for dark matter 
annihilating to 2e + 2e _ are illustrated in fig.[5]for two dif¬ 
ferent mass ratios m^/m x . The strong constraints from 
AMS lead us to models with negligible direct annihila¬ 
tions into electrons. 


3.5. Antiproton ratio 

The current measurements of the cosmic ray antipro¬ 
ton flux agree well with the expected astrophysical back¬ 
grounds; therefore, they are able to significantly con¬ 
strain the cross section for dark matter annihilation to 
hadronic final states, modulo the uncertainties from cos¬ 
mic ray propagation. Following [93], we use antiproton 
ratio data from the BESS [53] [55], CAPRICE [SB] and 
PAMELA [57] experiments. In order to minimize the 
systematic uncertainties from different experiments, the 
antiproton ratio data are employed here, rather than the 
antiproton flux data. Parameters governing cosmic ray 
propagation are constrained by fitting to the observed ra¬ 
tio of boron to carbon, but this still leaves some freedom 
to vary parameters, as well as the solar modulation; this 
leads to a range of estimates for the astrophysical back¬ 
ground of proton and antiproton fluxes, that are consis¬ 
tent with the observed fluxes from the above data sets. 
Once the propagation model, the solar modulation and 
the astrophysical background are fixed, we can perform 
a x 2 fit °f the p/p data by adding the contribution from 
the dark matter annihilation. 


1 The background may include contributions from primary and 
secondary electrons, secondary positrons, pulsars, or heavy dark 
matter annihilation/decay. All of them could contribute to the 
smooth spectra observed by AMS-02. 

2 The total magnetic field is composed of a regular one and a 

turbulent one. We normalize the total value at the Sun to Bq = 

15/tG. 
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Figure 4: Left column: distribution of versus m x for single-mediator models from fitting to CCW, Fermi and Daylan et al. 
data (from top to bottom), with f e = 0 and floating branching fractions to /r, r, q, b. Right column: corresponding distributions 
Of UJ-rJqJb- 


Two benchmark propagation models are used here 
to cover the systematic uncertainties of the cosmic ray 
propagation: KRA and THN, in the notation of ref. 
[92] , where the main difference of the two models is 
the halo height, with z t = 4kpc for KRA, and Zt = 
0.5 kpc for THN. The other parameters are the slope of 
the diffusion coefficient 5 = 0.5, nuclei spectral index 
7 = 2.35, the normalization of the diffusion coefficient 
Dq = 2.68 x 10 28 cm 2 s -1 , and the solar modulation po¬ 
tential >!> = 0.95 GV for the KRA model. For THN these 
take the same values except for Dq = 0.32 x 10 28 cm 2 s _1 . 
Neither of these models consider convection. Further de¬ 
tails can be found in |TKT ■ 

The resulting antiproton constraints on DM annihilat¬ 


ing to bbbb with two different ratios of mediator mass 
to DM mass are shown in fig. [6] The dependence on 
propagation model is the largest uncertainty, with THN 
(dashed lines) giving a factor of 10 weaker constraints 
than KRA (solid lines), while the dependence on the me¬ 
diator mass is relatively weak (red versus blue curves). 
We adopt the more conservative bounds from the THN 
model in our analysis. It will turn out that in our pre¬ 
ferred models with a subdominant branching fraction to 
b quarks, the antiproton limit is not significantly con¬ 
straining (see fig. [9] below). 

Better measurements of heavy nuclei abundances from 
AMS-02, in particular the boron/carbon ratio, may im¬ 
prove our understanding of cosmic ray propagation mod- 
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Figure 5: AMS exclusion curves for X + X > 20 —> 4e. The 
exclusion limit for two values of m$/m x = 0 and 0.9 are 
shown here, to illustrate the (relatively weak) dependence on 
the mediator mass. 



Figure 6: Antiproton exclusion curves for x+X —> 20 — > bbbb. 
The 3cr exclusion limit for two propagation models, KRA and 
THN, and for two values of m^,/m x = 0.5 and m < j,/m x = 0.8 
are shown here respectively. 

els, hopefully leading to a decrease in the systematic 
uncertainty from the propagation models. Moreover a 
much more precise measurement of antiproton data itself 
is expected to improve this limit by a factor of two or 
three. Considering both effects, we anticipate that com¬ 
ing AMS-02 data may shed further light on which DM 
models can consistently explain the GC excess. 

4. FITS WITH A SINGLE MEDIATOR 

Our initial motivation was to reexamine the viability 
of kinetically mixed Z' models for the GC excess in light 


Data set 

m x 

m<t, 

U 

fr 

u 

h 

(av) 

Xmin 

DOF 

CCW 

46 

12.3 

0.82 

0 

0.02 

0.16 

1.1 

18.8 

24 

Fermi 

130 

114.5 

0.80 

0 

0 

0.20 

2.8 

6.4 

20 

Ref. g] 

14.6 

4.0 

0.49 

0.01 

0.06 

0.44 

0.7 

22.9 

25 


Table I: Best fit parameters for single-mediator model fits to 
the three data sets. Masses are in GeV units, (av) in units 
of 1CT 25 cm 3 /s. “DOF” is the number of data points in each 
set. 


Data set 

m x 

rri0 

Xmin 

Data set 

m x 


(av) 

Xmin 

CCW (n) 

8.5 

i 

37 

CCW (b) 

75 

37 

0.6 

23.6 

Fermi (/r) 

14 

i 

42 

Fermi ( b) 

153 

102 

1.5 

20.5 

Ref. JH (n) 

7.5 

i 

51 

Ref. g] ( b ) 

51 

25 

0.4 

27.5 


Table II: Best fit parameters for single-mediator model fits to 
the three data sets, allowing for 100% branching to p final 
states (left) or b (right). Masses are in GeV units, (av) in 
units of 10~ 25 cm 3 /s. 


of the new data sets from CCW and Fermi. Although a 
good fit to the GC excess can be obtained, these models 
necessarily have a significant branching fraction into elec¬ 
trons. We found that the combined fit to the GC excess 
and AMS electron data was quite poor. 

This suggests doing a model-blind search for differ¬ 
ent combinations of final state fermions that could give 
a good fit to all data. We performed this in the 8- 
parameter space of models characterized by the DM and 
mediator masses m x , m^, and branching fractions _/) for 
annihilation into final states i = e, fj,, r, q , c, 6, where q 
denotes light quarks, whose spectra are all provided by 
PPPC [86]. These fractions are subject to the constraint 
yT fj = 1, where we ignore invisible channels since this 
degree of freedom would be degenerate with the overall 
cross section (av), which makes the 8th parameter. 

The result of this search is that f e must be negligible 
as a result of the AMS constraint. Therefore for more 
refined searches we set f e = 0 and remove the AMS con¬ 
tribution from the total x 2 - We find that the CCW and 
Fermi GC spectra further disfavor any significant contri¬ 
bution from r, q or c, preferring an admixture of muons 
and b quarks. We sample models with x 2 near the mini¬ 
mum value using a Markov Chain Monte Carlo (MCMC) 
search as well as Multinest [92]; this prevents getting 
stuck in spurious local minima of the x 2 function. In this 
way we find a large spread in allowed masses m x and m q > 
as shown in fig. [7] The best fit values of parameters for 
the three data sets are listed in table [J 

However the nominal best-fit values should not be 
given too much importance, because there are quasi¬ 
degeneracies in the x 2 functions that allow for good fits 
to the data over a range of parameters for the Fermi 
and CCW data sets. These are illustrated in fig. [7] We 
see that rn x can vary in the range 40-100 GeV keeping 
Ay 2 < 2 for the CCW fits. Similarly the Fermi data 
are compatible with m x in the range 80-180 GeV. Both 
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Figure 7: Top row and bottom left: best-fit values of m# (in GeV), / M and fb as a function of m x for data from CCW (top 
left), Fermi (top right) and Daylan et al. (bottom left). Bottom right: distributions of cross sections versus m x (1-2<t allowed 
regions) for the three data sets, compared to estimated Fermi limits from yy annihilation to p/r + bb in dwarf satellites (see 
text). 


data sets are compatible with annihilations into /i/i and 
bb , with the former being dominant. The corresponding 
results for the Daylan et al. spectrum are qualitatively 
different, both for the lower range of m x and the final 
states, which prefer /i, q or b depending upon m x . A 
range of larger masses m x ~ 55 — 70 GeV is however still 
reasonable with Ay 2 ~ 6. This provides overlap with the 
preferred regions of the other data sets. 

The required values of the annihilation cross section 
(av) needed to fit the magnitude of the excess are illus¬ 
trated in fig. [7] (bottom right). The l-2cr favored regions 
in the m x -(av) plane are somewhat disjoint between the 
three data sets. We show for comparison estimated limits 
on annihilation to 80% /ifj, + 20% bb from Fermi observa¬ 
tions of satellite galaxies of the Milky Way [100] , with an 
interpolation to the case of interest for us of / M ~ 0.8, 
fb ~ 0.2 using (av) -1 = f fJ .(a tJ .v)~ 1 +f b (a b v)~ 1 . Here the 
it’s refer to the limiting values. This is a rough estimate 
of the expected limit on our best-fit models to the CCW 
and Fermi GC excess data, which is seen to be still con¬ 
sistent with most of the perferred parameter space. Also 
plotted is our estimate of the limit on (a b v) from Fermi’s 
Pass 8 data m, which is lower than that of ref. m 


by a factor of 5. These results indicate some tension be¬ 
tween the dwarf limits and the GC excess signal. This 
can be relieved somewhat, still consistently with our fits, 
by reducing the fraction f b so that the signal is more 
/u.-like. Fig. [7] (bottom right) shows that the constraints 
into pure /.1 final states are significantly weaker. 

Since the allowed ranges for m x include values greater 
than mw and my, we have also performed Monte Carlo 
searches of the enlarged parameter space allowing for 
branching into W and Z. However we find that the best 
fit models have negligible branching into these states, and 
so we henceforth disregard them. 

One may wonder whether the preference for an ad¬ 
mixture of /i and b final states is statistically significant 
compared to 100% branching to either /j. or b. We have 
scanned over m x and for these cases to compare with 
the previous results for admixtures of final states. The 
best-fit parameters and minimum values of y 2 are given 
in table [TTJ to be compared to the results of table [T] An¬ 
nihilations to /i only give poor fits to all data sets, while 
those to b do better. Nevertheless, the improvement from 
allowing an admixture of /it and b is statistically signifi¬ 
cant. The q and r channels play a negligible role in the 



























fits to multiple final states, so these can be considered 
as approximately the same as having only mixed g and b 
states. The improvements in y 2 with two channels versus 
b only is 4 — 5 for the CCW and ref. gj data sets, and 14 
for the Fermi data. In all cases <5y 2 ^ 1 for the inclusion 
of one extra parameter. 

Physically, the preference for two final states arises be¬ 
cause the spectrum from any one of them falls with en¬ 
ergy as a power law times an exponential cutoff. Espe¬ 
cially for the Fermi data (see fig. [I]), the observed spec¬ 
trum does not have this form but instead remains high 
to large energies, just before abruptly cutting off. The 
superposition of final state spectra from g and b is able 
to reproduce this shape much better than that from bs 
alone. 


to fig. 0 (bottom right). Here we also show constraints 
on annhilations from Fermi dwarf observations (see dis¬ 
cussion of fig. [7]), and from our analysis of antiproton 
constraints, for the representative case of 20 % branching 
fraction to b quarks, except for the latest Pass 8 con¬ 
straint, where we use 10% to b quarks. There is strong 
tension between this limit and the preferred region for 
the GC excess of the Fermi data, while some compatibil¬ 
ity with the CCW data remains. The antiproton limit 
on the cross section is relaxed by a factor of 5 compared 
to the limit for annihilations purely into bs shown in fig. 
[ 6 ] The cross sections needed for the GeV excess are still 
comfortably below this limit (and more so in the single¬ 
mediator model where the target value of the cross sec¬ 
tion is lower). 


5. FITS WITH TWO MEDIATORS; 5.1. Sommerfeld enhanced annihilation 

SOMMERFELD ENHANCEMENT 


The previous analysis shows that annihilation to a mix¬ 
ture of b and g final states gives the best fits to Fermi 
and CCW data. From a model-building point of view, it 
may seem ad-hoc to couple the mediator to the standard 
model fermions in this way. If the mediator is a singlet 
that mixes with the Higgs, the couplings are proportional 
to the masses of the fermions. With two such mediators 
however, the desired mixture of final states could be a 
natural consequence of the masses m ^ and m^, 2 , only re¬ 
quiring that TO^q < 2 m T and 2 mb < m^, 2 < to x , so that 
decays into disfavored channels are kinematically forbid¬ 
den, or suppressed by small Yukawa couplings. Moti¬ 
vated by this theoretical consideration and the previous 
results, we thus reconsider the data in the model with 
five parameters 

{m x , m^, rrifc, (°" y )} 

This is one fewer parameter than in the single-mediator 
models, where we had two additional final states. 

Interestingly, the fits to Fermi and CCW data, in which 
m is free to vary, are consistent with values that are 
not far from the theoretical threshold 2 m T . In fig. [ 8 ] we 
plot best-fit parameters that are at least local minima of 
X 2 , demonstrating this assertion, y 2 is sufficiently flat 
as a function of m ^ that the models remain good fits 
even when m < f, 1 is restricted to stay below 2 m T , as is 
also shown in the figure. Monte Carlo searches of the 
parameter space reveal slightly better fits (<5y 2 ~ 0.1 — 
0 . 2 ) with higher ro^, but these do not contradict the 
goodness of fit of the models with m ( j 11 < 2 m T . 

Using MCMC to explore parameters in the model with 
fixed at 3 GeV, we find that the branching frac¬ 
tion to muons, f fl , is concentrated near 0.9 for fits to 
Fermi, while having a broader distribution in CCW. This 
is shown in fig. [9j left). The best-fit values of the annihi¬ 
lation cross section (crv) are shifted upwards by a factor 
of a few in these models relative to the single-mediator 
model, as can be seen from fig. [fright) in comparison 


The fact that these cross sections are significantly 
higher than the nominal value (crv )0 = 3 x 10” 26 cm 3 /s, 
needed for approximately the right relic density, is po¬ 
tentially a cause for concern. Such large cross sections 
would significantly suppress the density of y, in contra¬ 
diction to our assumptions. However the presence of the 
light mediator <f>± provides a possibility for resolving this 
problem, due to the low velocities of DM in the galaxy 
relative to the early universe, and the resulting nonper- 
turbative Sommerfeld enhancement of the cross section 
by multiple exchanges of <f>i. 

Let us suppose that (f>i couples to y with strength g\ 
and define op = <? 2 /47t. (Later we will see that both 
scalar g\ and pseudoscalar <71 5 couplings are needed to 
get s-wave annihilation, but that gi$ <C g±.) The Som¬ 
merfeld enhancement S is controlled by the two small 
parameters [ 102 * 


rrifa _ _ v 

€(b — 5 — 

aim x 


(7) 


A good approximation to S is given by the expression 

flMHIIMj 


S=%) - 8i ° llX ( 8 ) 

V e v J cosh X — cos yj(2ir/e ( j ) ) — X 2 

where = (w/12)e ( f, and X = Cy/e^. (The cosine be¬ 
comes cosh if the square root becomes imaginary.) 

To quantify the magnitude of Sommerfeld enhance¬ 
ment needed, and the corresponding coupling strengths 
a \, for each model in our MCMC chains we estimate the 
needed value of S as the actual cross section ( av) divided 
by the nominal relic density value (crv)o, assuming that 
v = 1CV 3 . Typically there are several coupling strengths 
that can yield the required value of S due to resonances, 
illustrated in fig. [To|left). Scanning over an, we find the 
minimum value of the coupling needed to get the required 
enhancement. The distributions of S versus ay for fitting 
the Fermi and CCW data are shown in fig. [To]( right). It 
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Figure 8: Left: best fits as a function of m x to the CCW data, as in fig. [7] Solid lines have m < / >1 freely varying, while dotted 
ones restrict m t j >1 = 3 GeV. Right: same for Fermi data. 



Figure 9: Left: distribution of branching fraction to muons for two-mediator model, with m t j >1 fixed at 3 GeV. Right: correspond¬ 
ing distributions of (av) versus m x , with Fermi dwarf galaxy constraints overlaid, as well as our constraint from antiprotons 
and that from the CMB. 


is noteworthy that we can get the right amount of en¬ 
hancement by invoking reasonably small values of the 
couplings, with 0.01 < aq < 0.06. 


6. PARTICLE PHYSICS MODELS 

To further explore the implications of dark matter 
with two scalar mediators, such as from the thermal relic 
density requirement and direct detection constraints, we 
need to more fully specify the interactions. We will allow 
for parity-breaking in the dark sector, with both scalar 
and pseudoscalar couplings 

2 

Ant = ^2 xMgi + igt, 5is)x (9) 

i -1 


the fact that the cross section for \X ~t 4>i(t>j is p -wave 
suppressed unless both types of couplings are present. 
For simplicity we take both mediators to be real fields. 

However, both scalars are required to get VEVs in or¬ 
der to mix with the SM Higgs and hence couple to SM 
fermions, so in fact eq. § would lead to a complex m x 
once the VEVs are taken into account, and we would 
have to redefine the fields to make it real. Instead we 
will regard ([9]) as describing the low-energy effective the¬ 
ory where this has already been carried out. 

The other important parameters of the model are 
the mixing angles 0; between the mediators and the 
SM Higgs. These are determined by cross-couplings 
\\i4i~\H\ 2 and VEVs, 


20 ; = 




mt - m% 


( 10 ) 


This choice is meaningful if the DM mass m x is restricted 
to be real; otherwise </; 5 can be rotated away by a chiral 
redefinition \ > e l6>75 x. The ansatz ([ 9 ]) is motivated by 


in the limit of 0; -C 1, where v = y/2{H) = 246 GeV. For 
our purposes, knowledge of the individual parameters A; 
and (<; t>i ) is not necessary, and we will work directly with 
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Figure 10: Left: example of Sommerfeld enhancement factor as a function of coupling ai = gr 2 /47r of light mediator 0i to dark 
matter, with m x = 100 GeV and m^ >1 = 3 GeV. Typical required value of S = 30 is shown by the horizontal line, for reference. 
Right: scatter plot of required values of S versus minimum value of ai corresponding to sample of models in fig. [9] 



Figure 11: 2<r-allowed values of mediator couplings to DM for fits to CCW (left) and Fermi (right) data sets. Regions differ 
depending upon the choice a o ne free combination of couplings which we take to be F (eq. (141, which we illustrate for F = 0 
and F = 0.9F max (see eq. (C4i). 


the 9i, that lead to couplings of the mediators to SM 
fermions of the form 

( 11 ) 

where yj = rrij jv is the SM Yukawa coupling. 

For simplicity we have neglected the mediator mixing 
mass m\ 2 0i02 and the cross-coupling Ai20i0l/4. They 
get generated at one loop by a virtual Higgs or y, at the 
level A 12 ~ (AiA 2 — ffiff!)/(167r 2 ). Cubic terms 0i0| 
and 020 2 can be forbidden by the discrete symmetry 
(j>i —» —0j, x e l7r75 / 2 x as long as x gets all of its 
mass from the singlet VEVs (in our case primarily (02)). 
This symmetry is spontaneously broken by the VEVs of 
0i so these terms get generated at one loop. They allow 
for the decay 02 —> 0101 which we have ignored. Taking 
it into account will require some small shift in the values 
of gi/g 2 needed to get the desired ratio of muons to b 
quarks for the GC excess signal, but since the data are 


not yet good enough to determine these couplings pre¬ 
cisely, we do not expect the 02 —> 0i0i decay channel to 
change our results significantly. 

To see how much fine-tuning is required by our ne¬ 
glect of 0 i -02 mass mixing, the most important term to 
consider is the y-loop diagram connecting 0i to 02 , and 
leading to m 2 2 ~ 5i<72TO 2 /(167t). If m\ and m\ denote 
the diagonal terms in the mass matrix, the lightest eigen¬ 
value is given by 


Taking m 2 ~ m x ~ 100 GeV, the second term is of order 
— 1 GeV 2 , similar in size to m 2 (we have assumed mi = 3 
GeV for our benchmark models). Therefore no extraor¬ 
dinary fine tuning seems to be required in our model to 
keep TOi small (of course we ignored here the Planck- 
scale hierarchy problem and took the heaviest threshold 
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in the hidden sector for the estimate of the loop contri¬ 
bution). An accidental cancellation between the terms in 
(12) could even help to explain the smallness of of m\. 


7. RELIC DENSITY 

The relic density of x is determined by \X 
summed over all possible final states. At kinematic 
threshold we find that 


(TV re \ 


E 

2 = 1,2 


4- 


8-7T (m| - to?/2) 2 

F 2 


167 t (m 2 — m|/4) 


(13) 


where for simplicity we approximated mi = 0 in the sec¬ 
ond line, and defined 


F = (9i92,5 + 32^1,5) + (9192,5 - 9291,5) 


9 

4m. 2 


(14) 


The couplings are constrained by the fits to the GC 
excess indicating that the number of muons to b quarks 
produced in the annihilations is given by 


N(g) = Njfa) (9i9i,s) 2 + F 2 /A 
N(b) N(<h) ~ (9292, 5 ) 2 + F 2 /4 


(15) 


For simplicity we here om itted the dependence upon 
m 2/ m x implied by eq. (13), but the more exact expres¬ 
sion is given in appendix [C| (eq. ©)• 

By demanding that (crv) matches the canonical cross 
section (av)o = 3 x 10 -26 cm 3 /s for approximately 
achieving the observed relic density, and also that (C2) 


gives the required ratio of g/b for a given model real¬ 
ization, we get two constraints on the couplings. The 
requirement of sufficient Sommerfeld enhancement gives 
a third, fixing the magnitude of g\. We can take the 
value of F in (14) as a free parameter that can be varied 


to explore the range of possible solutions for <7^77^5 char¬ 
acterizing viable models. Details of the algebraic solution 
for the couplings are given in appendix [C| 

We have carried out the above procedure for the mod¬ 
els in our Monte Carlo searches that give good fits to 
the GC GeV excess, to find the ranges of allowed values 
for the couplings. Scatter plots of <71,5, g 2 , 92,5 versus 
g 1 are shown in fig. m for fits to the Fermi and CCW 
data sets, within the 2cr-allowed regions. Here we have 
assumed solutions for gi that give the smallest values 
of the couplings (since a quadratic equation must be 
solved leading to a second branch of solutions with larger 
values). Interestingly, the parity-conserving couplings 
turn out to be the largest, with g\ ~ g 2 ~ (0.5 — 0.8), 
while the parity-violating couplings are suppressed, with 
91,5 ~ (0.02 — 0.1) and \g 2 ,5 | being smaller. The ex¬ 
act range of <72,5 (which we take to be negative, see ap¬ 
pendix is the least constrained of all the couplings, 


with <72,5 = 0 always being a possibility (corresponding 
to F = F max ), since we have freedom to impose this as 
an extra condition while satisfying the remaining physi¬ 
cal constraints. But in no case can \g 2 ,s\ be large while 
maintaining that g final states dominate over b for the 
GeV excess spectral shape. 


8. DIRECT DETECTION 


The cross section for scattering on nucleons (with mass 
m p ) is dominated by exchange of the light (j>i mediator, 


(91 0 i ym p ) 2 


(16) 


in the limit that m x m p . Direct detection limits put 
an upper bound on <71 9± that depends upon m x as shown 


in fig. 12 for fixed = 3 GeV. Here y is the Higgs- 
nucleon coupling, which we take to be y = 1.3 x 10~ 3 
(see ref. m for a recent review). For m x ~ 100 
GeV as suggested by our fits to the GC excess, this 
gives <7i0i < 2 x 10 -5 . On the other hand the Som¬ 
merfeld enhancement determined in sec tion ED demands 
that g ~ 0.6, so 9\ < 3 x 10 -5 . Eq. (10) then implies 
Ai(</>i) < 4 x 10 -3 GeV. We can eliminate the singlet 
VEV using its relation to the <f>i mass and self-coupling 
Ai by m r j >1 = \/2Xi((j)i). Taking m ^ = 3 GeV, we see 
that the direct detection constraint implies a hierarchy 
between the scalar quartic couplings, 


Ai 


A 


1/2 


< 2 x 10 


-3 


(17) 


This is not technically unnatural since the coupling Ai 
only receives multiplicative renormalizations. 



Figure 12: Upper limit on DM coupling times mixing angle 
of the light mediator (j>\ from LUX hsb] direct search. 
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9. COSMOLOGICAL CONSTRAINTS 

We consider here the impact of cosmic microwave back¬ 
ground constraints, and those coming from scattering of 
dark matter on itself in this section. We will show that 
the models treated in this work are compatible with cur¬ 
rent data, but are close to the limits, with interesting 
potential to address puzzles in the small scale structure 
of galaxies predicted by noninteracting cold dark matter. 


9.1. Cosmic microwave background 

Dark matter annihilations into charged particles dur¬ 
ing the epoch of reionization can impact the temperature 
and polarization fluctuations of the cosmic microwave 
background (CMB) |107 ]- [TT2l . leading to contraints on 
the annihilation cross section that are particularly strong 
for low DM masses, scaling linearly with m x . At this late 
era, the DM velocity is already sufficiently low that the 
magnitude of Sommerfeld enhancement occuring in col¬ 
lisions at the galactic center will be operative also for 
the CMB. Therefore we can directly compare the cross 
sections needed for the GC excess to the CMB limits. 

The CMB limits are depend upon the efficiency f 0 g for 
given final state particles to deposit energy electromag¬ 
netic energy in the plasma. The latest upper limit on the 
annihilation cross section using Planck temperature and 
polarization data can be expressed as 

fes{av) < 4 x lCT 26 ( 1QQ m Q eV ) ™ 3 /s (18) 

which we infer from fig. 41 of ref. m- Ref. dDI gives 
/ e ff = 0.25 (0.33) for g ( b ) final states at rn x = 100 GeV, 
which we average to f e g = 0.27 for our benchmark two- 
mediator model with 80% n + 20% b final states, to give 
(av) < 1.5 x 10 _25 (m x /100 GeV)cm 3 /s. This is plot¬ 
ted on fig. [9] (short dash-dotted curve). It is a some¬ 
what weaker constraint that the latest (Pass 8) limit from 
Fermi dwarf galaxy observations, but still excludes the 
preferred regions for the GeV excess fits to Fermi data 
with two-mediator models, while remaining marginally 
compatible with the fits to the CCW data. 


9.2. Dark matter elastic self-interactions 

Self-interactions of dark matter can be significant in 
our two-mediator model, from t- and u-channel exchange 
of the lighter mediator. The viscosity cross section, rele¬ 
vant for effects of DM self-scattering on structure forma¬ 


tion, is 3 

cr.„ = f dCl (1 - cos 2 8 ) = 9l m 4 x (19) 

J dil 3t r 

Taking the typical values g\ ~ 0.6, m ( / )1 = 3 GeV, 
rn x = 100 GeV, indicated by our fits, this leads to a cross 
section per DM mass of order a v /m x ~ 10 _5 b/GeV, 
far below the bound ~ 0.5b/GeV from simulations of 
structure formation including DM self-interactions cm. 
However, there is freedom in our model to take as 
small as 2771^ = 0.2 GeV, if we adjust op to somewhat 
smaller values ~ 0.01 to compensate for the increased 
boost from Sommerfeld enhancement in the galactic cen¬ 
ter. 

The self-interaction cross section in such models can 
have nonperturbative enhancements in analogy to the 
Sommerfeld effect, which have been studied in detail in 
ref. [1 1151 1116j . For aq = 0.01, the latter reference finds 
that a v /m x is at the right level to address problems of 
collisionless cold dark matter for predicting the observed 
galactic small scale structure, if m x = in the 

region = 0.2 — 2 GeV, which overlaps with values 
needed for explaining the GC excess. These problems 
include the difficulty for collisionless cold dark matter 
to correctly predict abundances and maximum masses of 
dwarf satellite galaxies, as well as the cusp versus core 
issue for dwarf galaxy DM density profiles; see ref. m 
for a review. 


10. COLLIDER CONSTRAINTS 

We also consider the collider limits on the particle 
physics model described in section [6] In a model hav¬ 
ing additional scalars, which mix with the SM Higgs, one 
would expect possible limits resulting from the recent ob¬ 
servation and measured signals of the Higgs boson at the 
LHC. Although we find the limits posed by direct detec¬ 
tion and relic density results to be more constraining, we 
include here the collider considerations for completeness. 

For the following anaysis, we do not initially assume a 
small mixing angle approximation. We supplement the 
SM Higgs potential with two real, scalar, singlet fields. 
Both fields acquire vevs, thereby inducing mixing with 
the Higgs. The relevant parameters are then the masses 
of the three scalars, and the two mixing 

angles, 81 hi 82 k- Mixing between the field and mass eigen¬ 
states is given according to 


3 for scattering of identical particles, this is more appropriate than 
weighting by (1 — cos0) since it treats scattering by 180° as 
equivalent to forward scattering 
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c 0 2 h 



( 20 ) 


where (c, s) denote (cos, sin), respectively. 

The scalar couplings to standard model particles are 
then simply the SM Higgs couplings, scaled by cg lh cg 2h , 
sg lh , cg lh sg 2h , for h, (f> i, and cf> 2 respectively. We take h 
to be the recently discovered Higgs boson, setting rrih = 
125.6 GeV, and take </>i to be the lighter scalar, such 
that m ( j) 1 < rrifo. We determine the allowed values of 
the two mixing angles, using the publicly available code 
HiggsBounds 4.2.0 pT£HHTl and HiggsSignals 1.3.0 H251 
1123) . We fix the values of the masses, such that rrih agrees 
with the discovered Higgs, as discussed above, and the 
values of rn ( f >i are those preferred by fits to the GC excess, 


m 0| =3 GeV, =115 GeV 



sin# lh 




in oj (GeV) 


Figure 14: Upper limit on </>i-Higgs cross coupling vs. m^ i; 
from CMS constraint on the total Higgs width, at 95% c.l. 
The shaded region is excluded. 


?7i0 1 = 3 GeV and m^ 2 = 115 GeV. 

The result is shown in fig. |13[ The upper figure shows 
the regions excluded by LEP and hadronic (Tevatron 
and LHC) exclusion Higgs searches, obtained by Hig¬ 
gsBounds, while the lower one shows the regions pre¬ 
ferred by compatibility with the observed Higgs signal 
strengths. The set of experimental results and specific 
Higgs channels that we use for the Higgs signal limit, 
can be found in refs |124Wl33j . For a complete list of the 
results used in the exclusion bounds analysis, see m 
and references therein. 

In the described analysis of Higgs signal limits, we only 
consider the effect of the mixing angles on scaling of the 
scalar couplings. In principle, new contributions to the 
width may further suppress the branching ratio, and sig¬ 
nal rate, in some channels. For simplicity, we present 
the results taking only SM decays of the scalars, as the 
collider limits turn out to be much less constraining than 
other limits considered here. 

For the mass hierarchy that is suggested by the GC 
considerations, one must also consider the additional con¬ 
tributions to the Higgs width, arising from decays to the 
light scalar, h —> 4>i(f>i. Returning to the small mix¬ 
ing approximation, the new contribution to the width is 
given by 


r 




\\v 2 

327rmft, 



( 21 ) 


Figure 13: Excluded regions of ( dih^ 2 h ) parameter space, 
under collider constraints. Top: blue (red) shaded region is 
excluded at 95% c.l. by LEP (LHC/Tevatron) exclusion lim¬ 
its. Bottom: preferred regions, compatible with LHC Higgs 
signal strength measurements. Red and blue curves corre¬ 
spond to 68% and 95% c.l. 


A recent CMS analysis [134) presents an upper limit 
on the Higgs total width, obtained from the ratio of mea¬ 
surements of off-shell and on-shell Higgs production and 
decay in the H —> ZZ channel. CMS finds T^ < 4.2 T® M , 
with T^ M = 4.15 MeV. We find the upper limit on the 
cross coupling shown in fig. |14| 
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11. CONCLUSIONS 

In this study we compared two new data sets, 4 charac¬ 
terizing the excess gamma rays from the galactic center, 
to predictions from models where dark matter annihi¬ 
lates into light mediators that subsequently decay into 
standard model particles. In contrast to models with 
direct annihilation into heavy quarks, these are more 
easily compatible with direct detection and cosmic ray 
antiproton constraints, as we have demonstrated. In a 
first approach, general admixtures of final state particles 
yielded good fits, revealing a preference for the data to 
be described by decays mostly to muons (electrons being 
excluded by AMS-02 data) and b quarks, with dark mat¬ 
ter masses in the range 40 — 140 GeV. Either Majorana 
or Dirac dark matter are viable possibilities. 

We then argued that this might be more naturally ac¬ 
complished with two scalar mediators <\>i mixing with 
the Higgs, such that <j>i decays primarily to /i + /U and 
(f >2 — > bb, as a consequence of the mediator masses. En¬ 
couragingly, fits to the data in this two-mediator model 
were consistent with masses in the desired ranges 2 tt< 
m ( f, 1 < 2to t , 2m,}, < m^, 2 < m x . Moreover quite reason¬ 
able perturbative values of the mediator couplings to \ 
give consistent results. To avoid p -wave suppression of 
the annihilation cross section we invoked parity violating 
couplings, which turn out to be somewhat smaller than 
the parity-conserving ones. 

An interesting feature of these models is that the best- 
fit cross sections for annihilation in the galactic center are 
several times larger than the value needed for achieving 
the right thermal relic density. We showed that this can 
be consistent due to Sommerfeld enhanced annihilation 
in the galaxy, due to multiple exchanges of the lighter of 
the two mediators, taking m^ 1 = 3 GeV as a benchmark 
model. It only requires that the scalar coupling of 
to % be of order <71 ~ 0.5. If the light mediator mass 
is somewhat lower, the cross section for elastic DM self¬ 
scattering can have the right magnitude for addressing 
the missing satellite and cusp-core problems from simula¬ 
tions of structure formation. Improved CMB constraints 
anticipated from Planck data may be in tension with the 
best fit regions of parameter space. 

LHC constraints on the couplings and mixing angles of 
the mediators are relatively weak compared to those from 
direct detection. The latter provides good prospects for 
independent confirmation of our model, requiring that 
the cross-coupling between <j> 1 and the Higgs boson be 
< 10~ 3 . No fine tuning in the technical sense is needed to 
satisfy this constraint, but neither is there any symmetry 
reason for the coupling to be small. 


The allowed parameter space in (<jv) versus m x shows 
some tension with the latest Fermi-LAT constraints on 
dark matter annihilation in dwarf spheroidals, especially 
for the fit to the Fermi GeV excess data. Planck con¬ 
straints on distortions to the CMB are in tension at a 
similar level. The fit to the CCW excess is also chal¬ 
lenged by these results, though to a lesser extent. More 
optimisitically, as we were completing this work, prelimi¬ 
nary evidence for a positive signal from the dwarf galaxy 
Reticulum II appeared |135j . at a level consistent with 
expectations from the GC excess. 
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E 1 [GeV] 

d$/d£ 7 dfi 

^stat 

^syst 

1.122 

1.587e-06 

1.036e-07 

8.225e-07 

1.413 

1.624e-06 

7.138e-08 

5.810e-07 

1.778 

1.483e-06 

5.330e-08 

3.240e-07 

2.239 

1.122e-06 

4.272e-08 

1.226e-07 

2.818 

7.298e-07 

3.655e-08 

5.857e-08 

3.548 

4.265e-07 

3.106e-08 

4.964e-08 

4.467 

2.475e-07 

2.074e-08 

3.511e-08 

5.623 

1.405e-07 

1.270e-08 

2.735e-08 

7.079 

7.662e-08 

8.267e-09 

1.874e-08 

8.913 

4.039e-08 

5.435e-09 

1.226e-08 

11.220 

2.272e-08 

3.688e-09 

7.959e-09 

14.125 

1.345e-08 

2.433e-09 

4.936e-09 

17.783 

7.828e-09 

1.566e-09 

3.016e-09 

22.387 

4.341e-09 

1.023e-09 

1.820e-09 

28.184 

2.503e-09 

6.953e-10 

1.115e-09 

35.481 

1.600e-09 

4.805e-10 

6.589e-10 

44.668 

1.029e-09 

3.146e-10 

4.090e-10 

56.234 

5.832e-10 

2.113e-10 

2.782e-10 

70.795 

2.753e-10 

1.355e-10 

1.556e-10 

89.125 

9.287e-ll 

7.851e-ll 

6.110e-ll 


Table III: Energy flux derived from Fermi collaboration’s pre¬ 
sentation [83] for the galactic center gamma ray excess, and 
statistical along with systematic errors as described in section 
|3.2| The flux intensity is obtained by averaging the observed 
total flux over the 15° x 15° square around the GC. Flux units 
are Ger'cm^s-'sr 1 . 


4 in particular, the Fermi collaboration spectrum determined 
within energy bands, as opposed to their ansatz using a power 
law with exponential cutoff, has not previously been analyzed 
with respect to dark matter models 


Appendix A: Fermi Spectrum 


We list the Fermi spectrum in table 
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Appendix B: Spectrum from decay to muons 


quarks resulting from the annihilations is 


For completeness, we present the photon spectrum 
from (j> —> y + in the rest frame of <f> [791187] , It includes 
photons from final state radiation and from radiative de¬ 
cays. Final state radiation gives the contribution 


d N FS R _ Otem 1 + (1 - 

da; 7r x 

where x = 2E 1 /m ( f,. 
decay is 


— 1 + ln 


1(1 — ^ 


( B1 ) 

The contribution from radiative 


dlV rac [ 
da; 


37 t x 


17 


191 2 23 _ 3 , 7^ 4 


--- -X + —rX - -X + -X 

2 2 12 3 4 


+ ( 3 + 7-x — 6a; 2 5 + 3x 6 — H 4 + 5a; In x ) In - 


+ ^3 + — 6a: 2 + 3a; 3 — Ml ~ x ) 

28 

——a; In a; + 5a;ln(l — a;) In a; + 5a;Li2(l — x) 

o 


(B2) 


where r = -C 1, and the range of x is (0,1) which 

does not depends on r since r is negligible. The sum of 
(Bll and (B2) gives us the total photon spectrum from 
annihilation to muons. 


Appendix C: Determination of couplings 


ax + cF 2 /2 
by + cF 2 /2 


(C2) 


Setting the cross sect ion equal to the relic density value 
(av)o and using (C2), we can solve for x and y, for a 
given value of F: 


x 


y 


1 f (Ho _ c p2 \ 

ayi + R- 1 2 ) 

1 ( (Ho c 2 \ 

b{TTR~ 2 F ) 


Since x, y > 0, this constrains 


(C3) 


R 

1 + R 
(C4) 

The minimum value of oq needed for sufficient Som- 
merfeld enhancement determines g\ = \Z4na 7; then 
Qi ,5 = \/x/g i and <72,5 = sfy/gii where we should al¬ 
low for both possible signs of y/x and yjy. Eliminating 
/i } 5 and f- 2,5 in equation |l4| for F results in a quadratic 
equation for s = g\jgi~- 


vHR’ 


F±^/F 2 -4^c^y(l-q 2 ) 
2 y/y(l + q) 


(C5) 


Here we show how the couplings of the mediators to 
dark matter are analytically determined from the various 
observational constraints. We parametrize the annihila¬ 
tion cross section as 


(<rv) = ax + by + cF 2 


(Cl) 


where x = (ffiffi.s) 2 , V = (5252,s) 2 , and o,6,c are the 
functions of rn x , m r p 2 in (13). The ratio R of muons to b 


where q = m^/4771 2 . The smallest value of |s| is found 
when the signs of the two terms in the numerator of (C5) 
are opposite. We can take F > 0 and the lower sign while 
keeping s > 0 if yjy < 0, corresponding to the choice of all 
couplings except for g 2 ,5 being positive. (Alternatively, 
we could take all couplings except for <7^5 positive and 
F < 0, but there is no physical difference.) One can 
explore the range of possible couplings by allowing F to 
vary between 0 and F max . 
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